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Abstract 

We study prediction-assimilation systems, which have become routine in meteorology and 
oceanography and are rapidly spreading to other areas of the geosciences and of continuum 



j^ ' physics. The long-term, nonlinear stability of such a system leads to the uniqueness of its 

sequentially estimated solutions and is required for the convergence of these solutions to the 
system's true, chaotic evolution. The key ideas of our approach are illustrated for a linearized 
Lorenz system. Stability of two nonlinear prediction-assimilation systems from dynamic me- 
teorology is studied next via the complete spectrum of their Lyapunov exponents; these two 
systems are governed by a large set of ordinary and of partial differential equations, respec- 
tively. The degree of data-induced stabilization is crucial for the performance of such a system. 
This degree, in turn, depends on two key ingredients: (i) the observational network, either 
fixed or data-adaptive; and (ii) the assimilation method. 



Physical systems — in nature, the laboratory or industry — can only be measured 
at a limited number of points in space and time. Estimating the state of a nonlinear 
dynamical system from partial and noisy observations is therefore crucial in applied 
physics and engineering [li[2j. In numerical weather and ocean prediction, this classi- 
cal estimation problem goes under the name of data assimilation [3l|4]; as data assimi- 
lation is spreading rapidly to other fields of the geosciences and of continuum physics, 
it is important to better grasp its fundamental theoretical aspects. In practice, a 
data assimilation algorithm has to provide the best-possible estimate of the evolving 
state of the system, using the observations available and the equations governing the 
system's time evolution [5]. In this paper, we examine the long-term stability of the 
set of modified equations that are referred to as the prediction- assimilation system, in the 
case in which the original physical system is fully nonlinear and chaotic. 

1 Introduction and motivation 

1.1 Background 

The complete solution of the filtering and prediction problem in sequential-estimation theory [HE] 
is given by the probability density function (PDF) of the unknown state, conditioned on the obser- 
vations. Given the correct initial PDF and assuming that the system noise and observational noise 
are white, normally distributed, mutually uncorrelated and known, the PDF's time evolution can 
be predicted by the Fokker-Planck equation [T] . In the case of a continuous stochastic dynamical 
system, with partial observations distributed at discrete times, an ideal data assimilation scheme 
would solve the Fokker-Planck equation for the time interval between observations and modify the 
PDF by using all observations when available. 

The fundamental difhculty of this approach relates to the high dimension of the state space, 
which makes it impossible in practice to obtain the initial PDF, let alone compute its time evolution. 
In the case of linear dynamics and of observations that are linearly related to the system's state 
variables, the PDF is fully characterized by its first and second moments, i.e. by the mean and 
covariance respectively. The optimal solution of the data assimilation problem in this linear setting 
is provided by the Kalman filter (KF) equations that describe the time evolution of both the mean 
and the covariance (TJ [21 [3l 131 [5l [6l [7] . The time-dependent error-covariance matrix depends, in 



this linear case, only on the observational error statistics that are part of the problem statement, 
and not on the actual observations [T] . 

Thus, in the case of linear dynamics, a linear observation operator, and observational and 
system noise that are both Gaussian, white in time and mutually uncorrelated, the KF equations 
give the optimal linear estimate of the state of the system by propagating the associated error 
covariances, along with the state estimates. In the nonlinear case, the situation is vastly more 
difficult and the PDF cannot be described by a finite set of parameters. A straightforward way of 
extending the linear results to the nonlinear case is given by the Extended Kalnian filter (EKF) 

[HEiiais]. 

In the EKF the tangent linear operator is used for predicting the approximate error statistics, 
while the state evolves according to the full, nonlinear equations. The computational cost of 
the EKF, though, is still prohibitive in many realistic circumstances. To alleviate this problem, 
a number of authors have studied reduced-rank approximations of the full EKF [S] UHl HH [12] 
that allow a reduction of its computational cost, while maintaining a satisfactory accuracy of the 
sequential estimates. A Monte Carlo approach, referred to as Ensemble Kalman filter, has also 
proven effective in reducing the computational cost associated with the full EKF [SI [131 [H] • 

Early theoretical work on the observability and stability of distributed-parameter systems (i.e., 
systems governed by coupled partial differential equations) was confined, by-and-large, to linear 
dynamics and to predetermined observations. In the case of linear, lumped-parameter systems 
(i.e., systems governed by coupled ordinary differential equations) and observations that are both 
linear and discrete in time, a sufficient condition for the KF solution to be stable is given by 
the observability of its dynamics [TJ [7]. Cohn and Dee [TS] have shown, in the linear, infinite- 
dimensional case of distributed-parameter systems that it is important to consider observability 
in the context of the discretized system, and that this observability implies stability of the data- 
assimilation problem. 

The concepts of observability and stability for nonlinear chaotic systems are closely related 
to other areas of dynamical system theory, namely controlling chaos and synchronization. In the 
control of chaos, a significant modification of the system's behavior is achieved by small variations in 
time of some parameter. Originally devised to stabilize unstable periodic orbits [16j . this approach 
has been generalized to force a given dynamical system to achieve other desiderable types of 
behavior, wheter stationary, periodic or chaotic [17] . In the present context, synchronization of 



chaotic systems means essentially using an adaptive coupling to have a "slaved" system track the 
motion of a driver or " master" system [T71 [THj . 

At the core of both chaos control and syncronization lies the stability problem. In the former, 
the time-dependent control has to be chosen so as to stabilize the motion. In the latter, the stability 
of the synchronized motion is a necessary condition for achieving such a motion. 

Interesting applications of both chaos control and synchronization to geophysical problems 
include the work of Tziperman et al. [19j on stabilizing an unstable periodic orbit in a fairly 
realistic El Nino model governed by a set of coupled, nonlinear partial differential equations, as 
well as that of Duane et al. [301 on meteorological teleconnections between the Atlantic and 
Pacific sectors of the Northern Hemisphere. Moreover, the relation between synchronization and 
data assimilation has been investigated by Duane et al. [31] and by Yang et al. [52] . 

1.2 The present approach 

In a chaotic system, initial errors grow within the system's unstable subspace. Trevisan and Uboldi 
[23| I24j considered fully nonlinear and possibly chaotic dynamics and proposed, in the context of 
meteorological data assimilation, to detect and eliminate the unstable components of the forecast 
error. They showed that those observations that help detect such instabilities maximize error 
reduction in the state estimates. Ghil i8j and associates (see references there) had already shown 
that, in meteorological and oceanographic data assimilation, the number of observations necessary 
to track an unstable flow is comparable to the number of the flow's dominant degrees of freedom, 
while Carrassi et al. [25 showed that one can improve on this estimate, since the requisite number 
of "tracking observations" is closely related to the number and magnitude of the system's positive 
Lyapunov exponents. 

In this paper, we examine the long-term stability of prediction-assimilation systems. This 
stability is essential for the performance of data assimilation methods and the convergence of 
their sequential estimates to the correct evolution of the underlying physical system. We de- 
velop a theoretical framework for the study of this long-term stability, and present a theorem 
that, under certain simplifying assumptions, provides rigorous conditions for the stability of the 
prediction-assimilation system. Within this framework, we present an approach that optimizes the 
convergence of the estimates to the correct solution, and apply it to two meteorological models of 
increasing complexity. 



The paper is organized as follows. In Sect. 2 we describe the formulation of data assimilation for 
nonlinear, chaotic dynamics, with particular emphasis on the proposed approach of Assimilation 
in the Unstable Subspace (AUS). Section 3 presents first the theorem, its proof and an illustrative 
numerical example; this illustration is followed by numerical results on the two nonlinear models, 
one governed by a large system of ordinary differential equations, the other by a system of coupled 
partial differential equations. Concluding remarks appear in Sect. 4. 

2 Data assimilation for chaotic dynamics 

We concentrate here on dynamical systems that are perfectly deterministic but chaotic; the role 
of explicit stochastic forcing, which may represent unresolved scales of motion, will be considered 
in future work. Without loss of generality, we write the system as a mapping from an arbitrary 
initial state at time io to a later time t: 

x(i)=X(x(io)), (1) 

where x is the n-dimensional state vector and M the nonlinear evolution operator. Given the 
initial state x(to) = xq, Eq. (1) will predict the state at future times t. However, due to the 
chaotic nature of the system, initial errors will amplify in time, thus setting a limit to the system's 
predictability. 

The tangent linear equations describing the evolution of infinitesimal perturbations (5x relative 
to an orbit of Eq. (1) can be written as: 

(5x(i) = M5x(io), (2) 

where M = M(x(io),t — io) is the linearized evolution operator associated with M, along the 
portion of trajectory between to and t. A chaotic system possesses one or more positive Lyapunov 
exponents, while their full spectrum characterizes the system's stability properties; these properties 
are crucial for the filtering, as well as for the prediction problem [U [26] . 

Suppose we seek an estimate of the state of this chaotic dynamical system from a set of noisy 
observations, given at discrete times tk >to,k^ {1, 2, ...}, 

y2.-H(xfc) + £^.; (3) 

here y^ denotes the p-dimensional observation vector, x^ the unknown true state, and eJJ is the 
observational error, all at time tk, while Ti is the (possibly nonlinear) observation operator. The 



observational error is assumed to be Gaussian with zero mean and known covariancc matrix R. 
We consider the underdetermined situation p < n; typically p << n in applications. 

To obtain an estimate of the state of the system, referred to in meteorological practice as the 
analysis x", one combines all available observations at tk with the background information, which 
consists of the forecast state at tk- This update is given by the analysis equation: 

x^ = (I-KfcW)x{+K,.y°., (4) 

where x^ indicates the forecast state, and Kfe is the gain matrix at time tk- We use here the 
unified notation [27] for meteorological and oceanographic data assimilation. In most sequential 
algorithms, the analysis equation has the form (4); such algorithms include the extended Kalman 
filter (EKF), as well as so-called optimal interpolation and other practical data assimilation schemes 
[UlUin]. Computing the optimally feasible K is at the heart of the sequential- estimation approach 
to filtering and prediction [U El [6] . 

The analysis state at time tk is obtained by applying the update (4) at this time to the forecast 
state x^ given by the nonlinear model evolution (1) of the analysis at the previous observation 
time tk-i'. 

x^ = (I-KfcH)X(x^i)+Kfcy°.. (5) 

The repetition of these analysis and forecast steps is referred to as the prediction-assimilation cycle. 

The effect of the observations can thus be interpreted as a forcing K(y° — HM{x°')), which acts 
on the free solution at the observation times tk', note that observations are typically not available 
at every time step of the discretized set of nonlinear partial differential equations O H] . Equation 
(5) governs the sequential estimation problem, i.e. the evolving estimate of the state of the system; 
y° — Ti.{A4{'x°')) here is the innovation vector. 

We consider now a perturbed trajectory that undergoes the same forecast and assimilation 
steps, with the same observations, as the reference trajectory of (5). The equation describing the 
linear evolution of perturbations dx-f'"'{tk) of this prediction-assimilation cycle is: 

SkI = (I - K,.H,.)Mfe_i<5x^_i, (6) 

where Hfe = H(x^,) is the Jacobian of TC at time tk and Mfc_i the linearized evolution operator 
associated with A4 between tk~i and tk. The term (I — KH) appearing in (6) reflects the effect 
of the forcing induced by the assimilation. This term modifies the stability properties of the 



perturbative dynamics relative to (5), i.e. its Lyapunov exponents, with respect to those of the 
free system (1). 

For the updates to drive the solution of Eq. (5) towards the correct solution of Eq. (1), 
the forecast-assimilation cycle (5) must be stabler than the pure-forecast system (1). Hence the 
Lyapunov exponents of Eq. (5) must be algebraically smaller than those of Eq. (1), which usually 
leads to its unstable subspace being lower-dimensional as well. Complete stabilization by the 
updating process, i.e. total absence of positive Lyapunov exponents, is sufficient for the uniqueness 
of the solution of Eq. (5), as well as necessary for the convergence of this solution to the true state 
of the system. Such complete stabilization will drive analysis errors to zero in the absence of 
observational and system noise, and to the lowest-possible values when noise or nonlinear effects 
are present. 

There are two means at our disposal in order to achieve this stabilization of a prediction- 
assimilation cycle: the design of the observational network, corresponding to the operator Ti, and 
that of the assimilation scheme, resulting in a certain gain matrix K; it is the product KH in 
Eq. (6) that provides the stabilizing effect of the forcing by the data. A. Trevisan and associates 
|23| , I24 1 [25] have proposed an efficient way to achieve this stabilization and improve the performance 
of the data assimilation method, by monitoring the unstable modes that amplify along a trajectory 
of the prediction-assimilation system. In their AUS approach, the basis of the subspace to which 
the analysis update is confined is given by the unstable directions of the system. 

The AUS gain matrix K, which differs from zero only on the unstable subspace at an update 
point, is given by: 

K = Er(HE)^[(HE)r(HE)^-|-R]"^ (7) 

Here E is the unitary matrix whose columns are the m unstable directions, while F is a symmetric, 
positive-definite matrix representing the forecast-error covariance in the subspace spanned by the 
columns of E; the index k is omitted in (7) for clarity. Since typically m << n, this feature of 
the method is clearly efficient in reducing the computational cost of the estimation process. A 
traditional, fixed network of observations can then be used to detect and reduce the forecast error 
projection along the unstable directions. An adaptive observational network, designed to measure 
primarily the unstable modes, will further enhance the efficiency of the assimilation. 

The unstable directions of a dynamical system of type (1) can be estimated by the breeding 
method [IHlllS]- In this procedure, the full nonlinear system is used to evolve small perturbations 



and, at fixed time intervals, their amplitude is scaled down to the initial value. The extension 
of the breeding technique to Eq. (5), referred to as Breeding on the Data Assimilation System 
(BDAS) [23], allows one to estimate the unstable directions of a prediction-assimilation system, 
subject to perturbations that obey Eq. (6). In the AUS assimilation, the unstable directions are 
used in the definition of the matrix K, cf. Eq. (7), and they can also be used to identify adaptive 
observations that are most beneficial for error reduction. 

3 Results 

3.1 Theoretical results 

We first provide a theoretical result that, under simplified circumstances, gives the mathematically 
rigorous condition for the observational forcing to stabilize the prediction-assimilation cycle (5). 
This result helps to clarify the theoretical underpinnings of AUS. Consider a chaotic flow, with a 
single positive Lyapunov exponent, and restrict the system's true evolution to an unstable fixed 
point, so that M in Eq. (6) is a constant matrix. The eigenvalues of this matrix are A^ — e^^^'^\ 
where A; are the Lyapunov exponents, and the eigenvectors of M are the Lyapunov vectors of the 
flow, while T = (tfc — tfc-i) is the assimilation interval. Alternatively, the result applies to the map 
associated with integer multiples of the period along an unstable periodic orbit of such a flow. 

Let the state of the system be estimated using a single noisy observation at each analy- 
sis time, assimilated by AUS along the single unstable direction ej,, so that: Kfe = CkBk and 



Cfe=72(Hefe) j^Ueky+a, 



v2 fXJr.. ^2 I ^2 



-1 
, 7^^ and a^ being the forecast error variance along e^ and the 



observation error variance, respectively. 

Theorem. Let the constant matrix M have a single eigenvalue A > 1 corresponding to a positive 
Lyapunov exponent, with e its associated eigenvector. Let H be a constant row vector, and 
the Kalman gain be approximated by K^ = cofe, where c is a constant scalar. The sequence 
{ej, : fc = 1, 2, ...} is defined by the recursion: 

Afeefe=Mffc_i, (8a) 

efcffc = (I-cefcH)Mffe_i; (8b) 

here the initial fo is an arbitrarily chosen unit column vector, while A^ and Qk are the normalization 
factors associated with e^ and ffc, respectively. Then, a sufficient condition for the solution of (5) 



to be stable is: 

cHe>l-A-i. (9) 

Remark. The amplitude c of the correction, and the observed (scalar) component He of the 
unstable vector, must thus be large enough to counteract the unstable growth. 

Proof. It can be shown that both e^ and f^ converge to e. The range of K is one-dimensional 
and the eigenvector e of M, with associated eigenvalue A, is also an eigenvector of (I — KH) M, 
with associated eigenvalue 

(I - ceH) Me = (1 - cHe) Ae = Ge. 

The stability condition is then obtained by setting 8 < 1, and this condition also guarantees 
uniqueness of the solution. D 

To illustrate the essence of the theorem we give here a simple numerical example in the context 
of the three- variable Lorenz model ^U\. When the canonical values of its parameters: a — 10, 
r = 28 and b = 8/3 are chosen, the system behaves chaotically. In the example that follows, the 
"true" state of the system that we want to estimate is the phase-space origin x ^ y = z = 0, which 
is an unstable fixed point for the system. Two assimilation experiments are performed; each one 
evolves according to the observationally forced, discrete dynamical system (5). At each analysis 
time a single noisy observation, the y- variable, is assimilated. 

The time step for integration has been set equal to 0.01 time units, while the assimilation is 
performed every 10 time steps, i.e., t = 0.1 time units; see also Miller et al. [31]. The matrix M 
represents the tangent linear operator evaluated at the origin and integrated for the time interval 
T. This matrix possesses exactly one eigenvalue larger than one, A — 3.26. The Kalman gain 
matrix, used to update the analysis, is K/j = ce^, with e^ obtained by the recursive use of Eq. 
(8). The sequence {e^} converges to the eigenvector e corresponding to the eigenvalue A, and the 
asymptotic value of He^ is He — 0.91. The stability condition (9) yields therefore c > cq = 0.762. 

Figure 1 shows the root-mean-square (RMS) analysis error as a function of time (in units of 
assimilation interval, i.e., r = 10 time steps) for the two experiments, with c = 0.76 and c = 0.77, 
respectively. The observational RMS error is 10^'^. Clearly the RMS grows exponentially in the 
former case and decays exponentially in the latter. 



3.2 Numerical results 

We now illustrate the stabilizing effect of observational forcing for different assimilation schemes 
and show that their long-term performance is related to the induced degree of stabilization. The 
full Lyapunov spectrum is used to compare the properties of a given free system (1) with those 
of the corresponding prediction-assimilation system (5). Observing system simulation experiments 
[21 m [5] are performed with numerical models of increasing complexity. 

The first chaotic model [32] has 40 scalar variables that represent the values of a meteorological 
field at equally spaced sites along a latitude circle. It can be derived from an "anti-Burgers" 
equation in the same way that the Fermi-Pasta-Ulam model was subsequently shown to be derived 
from the Korteweg-deVries equation |33j. The second is an atmospheric model that represents 
mid-latitude, large-scale flows, and is based on the quasi-geostrophic equations [9J in a periodic 
channel. The discretized model has 15 000 scalar variables and its details can be found in [51] , 
while the experimental setup for assimilating data is described in |25j . 

Figure 2 shows the spectrum of the 40 Lyapunov exponents of the first model (32] , using 200 
years of simulation. The free system possesses 13 positive exponents (A+), of which the leading one 
{^max = 0.336 day~^) corresponds to a doubling time of 2.06 days; its Kaplan- Yorke dimension is 
approximately 27.05. At each assimilation time, the analysis is performed by AUS with a single 
BDAS mode e^^"^^ in the specification of the gain matrix, Eq. (7). In this single-observation 
situation, the matrices E and T reduce to a column vector and to a scalar, respectively; the latter 
is estimated statistically using the innovations — which are scalars, in the present case — following 
the approach of ^2Sj . 

A single observation, adaptively located where the current BDAS mode attains its maximum 
value, is sufficient to stabilize the system, so that Xmax = —0.283 day^^, and to reduce the RMS 
analysis error to 1.4% of the system's natural variability, even when the assimilation interval 
T = tk+i — tk is as long as 3 hr. Smaller assimilation intervals of 2 hr and 1 hr lead to further 
stabilize the system {Xmax = —0.437 day~^ and —0.809 day~^, respectively), and thus to reduce 
the analysis error even more: the RMS error, normalized by the natural variability, becomes 0.011 
and 0.009, respectively. 

Figure 3 shows the first 100 Lyapunov exponents of the quasi-geostrophic model [M], using 
one year of simulated time. The free system possesses 24 positive exponents (A+), of which the 
leading one {Xmax = 0.310 day^^) corresponds to a doubling time of 2.2 days, and its Kaplan- Yorke 
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dimension is approximately 65.2. The three assimilation experiments all use a fixed network of 
noisy observations that cover just under one-third of the domain (20 out of 64 meridional lines 
of grid points); in two of the experiments, an additional observation is adaptively located at a 
single grid point in the otherwise unobserved portion of the domain. Its location coincides with 
the maximum of a single BDAS mode, as in the first model i32j. The model time step is At = 30 
min, while the assimilation interval r is 6 hr. 

In all three experiments, the fixed observations are assimilated by a least-square fit, according 
to the three-dimensional variational (SDVar) algorithm [SJ H] in wide operational use, while the 
adaptive observations are assimilated either by SDVar (3DVar-BDAS) or in the unstable subspace 
(AUS-BDAS) by using the current BDAS mode in Eq. (7). When fixed observations only are 
assimilated (SDVar), the number of positive exponents is reduced to three, with the leading expo- 
nent {Xmax = 0.088 day^^) corresponding to a doubling time of 7.9 days, while the Kaplan- Yorke 
dimension is reduced to 6.9 and the normalized RMS error to 0.S21. Adding a single adaptive ob- 
servation assimilated by SDVar (SDVar-BDAS) stabilizes the prediction-assimilation cycle further, 
the only positive exponent being slightly greater than zero {Xmax — 0.002 day^^, Kaplan- Yorke 
dimension 1.1), and reduces the normalized RMS analysis error to around 0.16. Finally, when the 
adaptive observation is assimilated in the unstable subspace, the system is completely stabilized 
and the RMS analysis error drops to only 0.058. 

4 Conclusion 

To estimate the efficiency of prediction-assimilation systems and observational networks, one of- 
ten estimates the error of a short-range forecast at points where fairly accurate observations are 
available; the obvious drawback of this approach is that errors tend to be smaller in systematically 
observed regions O [H [35] . The nonlinear stability analysis introduced here allows one to address 
these issues in a more rigorous way. The stability of the prediction-assimilation system guarantees 
the uniqueness of its solution and is required for the convergence of this solution to the true flow 
evolution; in turn, the degree of stabilization introduced by the data assimilation may be measured 
precisely by estimating the full Lyapunov spectrum of the forced system. 

Assimilation in the unstable subspace and the use of breeding on the prediction-assimilation 
system to estimate this subspace were applied here to a 40-variable [S2 and to a 15 000-variable 
model 25 simulating the mid-latitude atmospheric circulation; in both models, the proposed 
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methods led to complete stabilization of the sequential estimation process. The tools of dynamical 
systems theory can thus help design optimized assimilation algorithms (the specification of K) 
and observational networks (specifying Ti.) for high-dimensional and highly nonlinear systems. 
When, as usual in practical applications, only a limited number of measurements can be made, 
asymptotic error reduction may still be achieved through adaptive deployment of observations and 
a sophisticated assimilation scheme, designed to control the flow's instabilities. 

Data assimilation applications are possible in all situations where a dynamical constraint is 
important and only a limited amount of noisy observations can be taken. Such situations include 
robotics, fiow in porous media, plasma physics, as well as solids subject to thermal and mechanical 
stresses or to shocks [3Sl[3n]- The insights into the nonlinear dynamics of the prediction-assimilation 
cycle provided here should help design data assimilation methods for a large class of laboratory 
experiments [37], as well as for natural or industrial systems. 
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Figure 1: RMS analysis error as a function of time for two experiments with c = 0.76 and c — 0.77, 
below and above the stability threshold cq — 0.762; the time is given in multiples of the updating 
interval r = 0.1. 
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Figure 2: Effect of the observing system set-up on stability of the prediction-assimilation cycle 
for the Lorenz 40- variable model |_32i . Spectrum of the Lyapunov exponents for the free system 
(1) (diamonds) and for the AUS-BDAS forced system (5), given different assimilation intervals 
T = tk+i — tk'. 3 hr (squares), 2 hr (triangles) and 1 hr (circles). 
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Figure 3: Effect of the assimilation method on stabihty for the quasi-geostrophic model [51] . 
Spectrum of the first 100 (the first 10 in the inset) Lyapunov exponents for the free system (1) (di- 
amonds) and for the observationally forced systems (5): 3DVar (squares), 3DVar-BDAS (triangles) 
and AUS-BDAS (circles). 
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